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Abstract 

We study the modes and stability of non - isothermal coronal loop models with different intensity 
values of the equilibrium magnetic field. We use an energy principle obtained via non - equilibrium 
thermodynamic arguments. The principle is expressed in terms of Hermitian operators and allow 
to consider together the coupled system of equations: the balance of energy equation and the 
equation of motion. We determine modes characterized as long - wavelength disturbances that are 
present in inhomogeneous media. This character of the system introduces additional difficulties for 
the stability analysis because the inhomogeneous nature of the medium determines the structure of 
the disturbance, which is no longer sinusoidal. Moreover, another complication is that we obtain a 
continuous spectrum of stable modes in addition to the discrete one. We obtain a unique unstable 
mode with a characteristic time that is comparable with the characteristic life-time observed for 
loops. The feasibility of wave-based and flow-based models is examined. 

PACS numbers: Valid PACS appear here 
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I. INTRODUCTION 



A. Variational principles 



Stability is a crucial requirement for a model to produce realistic descriptions. Thus, dif- 
ferent stability analyzed of solar structures can be found in the literature, generally restricted 
to special types of perturbations and specific equilibrium models. These includes, models 
that consider adiab atic c onfiguration such as the ones analyzed via the classical criterion 
of Bernstein et al. 19581 ] or those that presuppose static equilibrium and analyze thermal 



stability. In the application of Bernstein's criterion, the adiabatic assumption implies that 
the energy balance equation is not required and thus dissipation is impossible. Also the 
assumption of static models is a strong, and often unjustified, restriction for open systems. 
Thus, a crucial question for any theoretical model is whether the much more common far- 
from-equilibrium states are stable, where the consideration of both thermal and mechanical 
coupled equations must be included. 

A more realistic analysis of the stability of config urations represented by non- conservative 



equations was presented by Lerche and Low 19811 ] . They proposed a Lagrangian principle 
in order to analyze quiescent prominences that can undergo thermal instabilities. However 
the non-self-adjoint character of the operators involved in the obtained principle makes the 
physical interpretation difficult. 

In this paper we apply an energy principle to analyze the stability of so 
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et al. 



ar coronal loops. 

e was obtained in a previ ous p a per (Paper I : Costa et al. |2004| ; see also Sicardi 



20041 ]. see also Sicardi et al. [1985 



1991 



1989a 



1989b|) using a general procedure 



of irreversible thermodynamics -based on firmly established thermodynamic laws- that can 
be understood as an extension of Bernstein's MHD principle to situations far from ther- 
modynamic equilibrium. This fact has the advantage that many known results obtained by 
simpler criteria can be re-examined by a direct comparison with our criterion, and that, as it 
is obtained via a thermodynamic approach, allows a straightforward physical interpretation. 
The principle associates stability with the sign of a quadratic form avoiding non-self-adjoint 
operators. Obtaining a self-adjoint operator is a requirement for our principle to hold. When 
this is accomplished the calculus is simplified. The self-adjoint character of an operator im- 
plies that the eigenvalues u 2 are real. Hence stability transitions always occur when u 2 
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crosses zero, rather than at particular points of the real axis where the real part of the 
eigenvalue is different from zero, i.e. Re(u) 7^ 0, leading to an efficient formulation to test 
stability. Thus, the symmetry considerations of the self-adjoint operators, the fact that 
there is a diagonal form associated with these operators, and that the Rayleigh-Ritz theo- 
rem states the existence of a minimum eigenvalue, are important reasons to try to maintain 
self-adjointness in the consideration of stability. 



1. Solar coronal loops 



The theoretical modeling and the interpretation of observations of coronal loop systems 
deal, among others, with the discussion whether the propagating observed disturbances in 
loops and post-flare loops are waves or plasma flow. 

Dynamical features of brightening coronal loops have been traditionally interpreted as 
field-aligned flow of matter generated by asymmetries in the energy input. Most classical 
theoretical models have difficulties in determining the physical conditions that make them 
compatible with observations. Both static loops and steady state models -for the two classes 
of temperatures models: hot (isothermal coronas with T pa ItfiK) and cool (gradually 
increasing temperatures up to T pa 10 5 K) - fail to provide a satisfactory explanati on fo r 



both t he em ission measure distributi on and the Dopp 



er sh ift observations (Jordan 19801 ]; 



Serio 19811 ] : Craig and McClymont 19781 ] ; Mariska 1984J). Thus, this suggests that in 
traditional model scenarios radiative losses cannot be compensated by therm al con d uction . 
Therefore, other he ating mechanisms must be assumed (Aschwanden et al. 1999] ; 2000]; 



Walsh and Galtier 



20001 ]). Also, theoretical time-dependent models of individual loops 



where the plasma evolves in response to a cyclical process of he ating and cooling of the flow 
have difficulties in fitting observations (Klimchuk and Mariska 1988] ). 

The assumption of propagating disturbances associated with slow magnetoacoustic waves 
in high Alfven speed media is also a field of investigati on. S everal wave-based models 
were d eveloped to explain observations (Nakariakov et al. 2000] ; Tsiklauri and Nakariakov 



2001] ). These authors suggest that -depending on the relative importance of dissipation 



by magnetic resistivity- upwardly propagating waves (of observed periods between 5-20 
min) that decay significantly in the vicinity of the loop apex could explain the rarity of 
observational detection of downwardly propagating waves. However, upwardly propagating 
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disturbances with non-decaying or even growing amplitudes were observed in coronal EIT 
plumes. Analytic models have shown that slow magnetoacoustic waves may be trapped and 
nonlinearly st eepen ed with height, providing a possible interpretation of the phenomenon 
(Ofman et al |1999| ). 

However, due to the intensity of the flaring, the plasma dynamic of flare loops is generally 
associated with flows rather than with waves. In fact, systematic intensity perturbations 
in post-flare loops can suggest that they are the result of evaporation-condensation cycles 
caused by the efficient heating of the flaring plasma from the chromosphere. Thus, chro- 
mospheric ev apora tion seems to be the main initial matter inflow source for flare loops. De 



Groof et al. [2004] analyzed an off-limb half loop structure from an EIT shutterless cam- 
paign and gave arguments to reject the slow magnetoacoustic description and to support the 
flowing/falling plasma one. Nevertheless, these authors admit that the wave theory cannot 
be excluded yet. 

Other authors have suggested that a comb ination of phenomena can be at the basis of a 
better interpretation. Alexander et al. 19981 ] examined 10 flares and concluded that plasma 
turbulence could be the source of the observed intensity changes rather than hydrodynamic 
flows such as chromospheric evaporation. They pointed out that it cannot be excluded that 
there is a degree of "gentle evaporation" occurring early in the event with associated hard 
X-ray emiss io n below t heir t hreshold of detection. A series of more recent papers (Tsiklauri 
et al. 20041 ] a; 2004j]b; 2004j]c) that combine theoretical and observational analysis showed 
that oscillations in white, radio and X-ray light curves observed during solar and stellar 
flares may be produced by slow standing magnetoacoustic modes of the loops. They found 
that a transient heat deposition at the loop bottom -or at the apex- leads to a posterior 
up-flow evaporation of material of the order of a few hundreds of km/s^ 1 . During the peak 
of the flare, the simulations showed that a combined action of heat input and conductive and 
radiate losses could yield an oscillatory pattern with typical amplitudes of up to a few tens 
of km/ s~ l . Then, a cooling phase of plasma draining with velocities of the order of hundreds 
of km/s^ 1 occurs. The numerical quasi-periodic oscillations in all the physical quantities, 
that resemble observational features, were interpreted as being produced by standing sound 
waves caused by impulsive and localized heat ing. 

In previous papers (Borgazzi and Costa [2004J; Costa and Stenborg [2004J) one of us 
developed a diagnostic observational method to describe loop intensity variations, both 
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in space and time, along coarse-grain loop structures. We find that none of the arguments 
leading to the determination of whether waves or flow models can better fit observations was 
conclusive. Some of our results suggested wave-based model interpretations i.e. the periodic 
behaviour of the disturbances observed, the almost constant speed of some brightening 
features and the fact that the estimated speeds were not higher that the sound speed in the 
coronal loops. Howe ver, a s we mentioned, t he pe riod behaviour can also be attributed to 
flows (Gomez et al. 1990 ]; De Groof et al. ^OoJ ]). Also, even when the calculated speeds 
were not greater than the sound ones, some of the velocity patterns were far from being 
constant and their values were comparable to the free-fall case. 

Another open question is the relation between the loop's coro nal dy namics and the phys- 
ical conditions on the chromospheric bases. Borgazzi and Costa 2004| found a longitude of 
chromospheric coherence that characterizes the behaviour of a whole loop-system of evolv- 
ing coronal-isolated filaments. This description is in accordance with limit-cycle models 
that require that the triggering mechanism of the dynamics is located at the bottom of the 
structure giving rise to the observed similar coronal conditions of the isolated filaments. An- 
other aspect that deserves attention is whether it is physically possible that the periodicity 
observed could be related to, or could be the consequence of propagating magnetoacoustic 
modes from the chromosphere that have suffered distortion due to the dispersing media. 

Other point that is under debate is the thermal structure of the loops . Loo p observations 
with TRACE (Transition Region and Coronal Explorer, Handy et al. [19991]) suggest that 
hot coronal loops are isothermal and more dense than the predictions of static loop models. 
Howe yer th is scenario is not conclusive and other interpretations are possible. Reale and 
Peres 2000| showed that bundles of thin strands, each one behaving as a static loop, with its 



characteristic thermal structure, convoluted with the TRACE temperature response could 
appear as a single almost isothermal loop. A wide range of configurations can be proposed to 
fit observations. The fact that images form a compound of complex integrated time-varying 
data that are not easy to resolve is at the basis of this difficulty. The loops under analysis 
are surrounded by other structures that usually intersect them along the line of sight and 
the change of the brightening of the loops is also affected by background emission. Thus, 
efforts are made to produce observational and theoretical results of coronal loop dynamics. 

The aim of this paper is to investigate whether the propagating observed disturbances 
in loops are waves or plasma flow and their thermal structure. Non-isothermal loops are 
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traditional candidates for Hopf instabilities with cycles of flow evaporating and condensing, 
thus the analysis of frequencies and mode structures can provide insight into a possible wave 
model interpretation of these types of configurations. We consider the stability analysis as 
the leading criterion to select possible theoretical wave models. The fact that a number of 
non-linear equilibria are possible due to the open character of the systems makes it necessary 
to consider both thermal and mechanical stability in a coupled way. 



II. THE STABILITY CRITERION 



The thermodynamics of irreversible processes is described in terms of phenomenological 
relations between conjugate pairs of thermodynamic variables: the flows and the forces that 
cause them. The linear thermodynamic approximation treats small deviations from the 
equilibrium state by including fluctuations in the neighborhood of this state. It describes 
the behaviour of the system around the equilibrium state or around a non-equilibrium 
stationary one that is linearly close to it. 

If the system is isolated, as is stated by the second law of thermodynamics, the en- 
tropy grows exponentially up to its maximum value. That the system is in an open-near- 
equilibrium state means that energy and matter is exchanged with the neighbors and the 
entropy of the system is not necessarily positive. Even when the entropy produced in the 
system's interior, due to irreversible processes, is never negative, a negative flow of entropy 
produced by the exchange of matter and energy can make the system remain indefinitely 
in a near-equilibrium state. These states are known as stationary states and a coherent 
dynamic of the system could last if sufficient negative entropy flow is provided to it. Thus, 
the criterion that states the stability of this stationary state gives insight into the dynamic 
structures that can be found in na ture. These stationary states are also known as detailed 
balanced. As Onsager pointed out 19311 ] . the balance consists of the compensation between 
the fluctuations and dissipation produced by the flows and forces that have a microscopic 
reversible character near the thermodynamic equilibrium. The empirical relations between 
flows and forces are linear and related by the so-called resistance matrix R that is symmetric 
and positive definite. Its symmetric character is guaranteed by the principle of microscopic 
reversibility and its positive definiteness by the proximity of the reference state to the ther- 
modynamic equilibrium, where the entropy has a maximum. 
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However, there is no continuity between linear and nonlinear thermodynamical processes. 
When the system is beyond the immediate neighborhood of the stationary state the non- 
linearities become visible. Instabilities that cause dynamic transitions in open systems 
are responsible for the qualitative difference between linear and nonlinear thermodynamics. 
Therefore, dynamic cooperative phenomena can only arise in nonlinear thermodynamics. 
Thus, nonlinear thermodynamics is related to the stability properties of non-equilibrium 
stationary states, where the linear relation between flows and forces can become state de- 
pendent (i.e. R is not necessarily a symmetric positive definite matrix), and the problem 
of having a thermodynamic theory to provide a general criterion for the stability of the 
system -which is not evident through the integration of the variational equations- becomes 
a fundamental point. Non-linear thermodynamics is the extension of the linear theory to 
situations far from thermodynamic equilibrium where the relaxation of the processes to a 
steady state of non-equilibrium (nonl inear sta te) is not assure d and requires a stability ana l- 
ysis (Glansdorff and Prigogine 197l| : Keizer 1976 ]; Graham 1978 ]; Lavenda 1987 . 1993 ]). 

In Paper I we showed how to obtain the variational principle from the equations that de- 
scribe the dynamics of the system of interest. The method consists of obtaining a Lyapunov 
function, also known as generalized potential, that represents the mathematical expression 
of the stability conditions. This function is determined by the analysis of the thermody- 
namic properties of the system linearized around a non-linear stationary state also called 
non-linear equilibrium state. The equations governing the dynamics are written as a system 
of two coupled equations: the balance energy equation and the equation of motion. Thus, 
the perturbation analysis around a stationary state is performed considering a variable state 
vector of four independent components: the three space component displacement and the 
temperature variation. Once the linearization is done, the Lyapunov function can be im- 
mediately obtained by inspection of the resulting expression written in a compact matrix 
form. Each of the matrices of the compact expression are linear operators (that could in- 
clude spatial derivatives) and have a clear physical interpretation that is given by its role in 
the equation. The matrix that multiplies the second time derivative of the perturbation is 
associated with the inertia of the system, the one that multiplies the first time derivative of 
the perturbation is associated with dissipation and the one that multiplies the perturbation 
itself is associated with potential forces over the system. The principle is subject to phys- 
ically reasonable requirements of hermiticity and antihermiticity over the matrices. For a 



7 



more detailed presentation see Paper I and the references presented there. 



A. The magnetohydrodynamic expression 

The specific model we analyze is taken to be composed of a magnetohydrodynamic ideal 
plasma (i.e. with infinite electrical conductivity cr 3> 1). The fundamental ideal magneto- 
hydrodynamic equations to be considered are as follows. The mass conservation equation, 

^ + v-(H = o (i) 

where p is the density of the plasma, v is the plasma velocity, and t the time. The perfect 
gas law or state equation, 

V = ~PT (2) 
m 

ks is the Boltzmann constant, p the pressure, T the temperature and m = m p the proton 
mass. For a fully ionized H plasma p = pn e m p ; the solar coronal abundances (H : He = 



10 : 1) correspond to a mo 



ecular weight p = 1.27; n e is the number density of electron 



particles (Aschwanden [2004f ] . The induction equation, 

— = Vx(vxB) (3) 

B is the magnetic field vector. The magnetic diflusivity was discarded. The equation of 
motion for the problem is: 

P Wt = ~ V(P) + ^(V x fi) x - P V0 (4) 

where g = — V0 is the gravity expression and j = l/47rV x B the current density. The 
energy balance equation takes the form: 

where 7 is the ratio of specific heats and L is the energy loss function: 

L = —V ■ F c — L r + H. (6) 
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F c is the heat flux due to particle conduction along the loop, L r is the net radiation flux. 
Neither the dominant heating mechanism of coronal loops nor the spatial distribution func- 
tion of the energy input is known. So, the heating function is the least known term in the 
energy balance equation. Thus, the usual situation is to try reasonable arbitrary mathe- 
matical functions which must fit the constraint imposed by the equilibrium conditions. As 
our model considers inhomogeneous temperature gradients and isothermality is usually as- 
sociated with footpoint heated loops more than with uniformly heated ones we discard the 
first case and tried the general expression H = hp + Hq. A time varying dependence of H 
was not considered for simplicity. However, it could be a requirement for modeling special 
events such as micro-flares or while considering magnetic reconnection phenomena. Eq. [5] 
expresses the fact that the gain in particle energy (internal plus kinetic) is due to heating 
sources, heat flow and radiation losses; ohmic dissipation j 2 /a and all other heating sources 
were considered as vanishing terms implying that the optically thin assumption holds. Then 



L r = n e riHQ(T); the temperature variation (Q(T) = xT a ) was taken from Priest [1982 ]. 
Also F c = —kVT and, as conduction across the magnetic field has been discarded, for a 
total ionized plasma F c = — k T^X7\\T. Finally equation [5] can be written as 



V • (k T* V,|T) - ^ X T a + ^P + H (7) 



7 — 1 Dt m 2 m 2 

where v is a constant value to be determined from the equilibrium conditions. 

The linearization procedure is performed by replacing p — p + p x , T = T + T 1; B — 
Bq + Bi and v = vq + d^/dt in the last equations, and assuming hydrostatic conditions 
for the equation of motion. Thus, vq = and V\ = d^/dt where £ is the perturbation 
around the equilibrium of the equation of motion (the stationary state), also dpo/dt = 
and dBo/dt = 0. Using the relation d/dt ~ iu in eq. [Hand eq. El the corresponding 
linearized equations (eq. M - H2]) are: 

Pi + V • (poO = (8) 



Pi = — (poTi-ToV-Cpofl) (9) 
m 
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B, = -V x (B x I) (10) 

p f=M V (T V-(p 3-Po7i)- 

_ J_[(v x Q) x £? + (V x S ) x Q] + V 0V • (pol) (11) 
4p 

or equivalently p |* - F£ + ^V(p Ti) = and 

.J PoTi - ( 7 - 1)T V • (Pol)] - ^Ti + Sf = (12) 
to (7 — Ij 



being 



A = -[ C V • (T *V||(*) + ^Vn(To)) - ^X«T?- ]1 



and 



m Kbtti 



1 810~ 10 

c = — — — Wm^K- 1 , Q = B x£. 
InA 



The term V- (po£) was discarded because it represents the total net flux of material through 
the magnetic tube. The two obtained equations are expressed in terms of the displacement 

— * 

and temperature perturbed variables £ and T\. * represents the location of the perturbed 
variables when performing the matrix product. 

Following Paper I the resulting energy principle is: 

+ J (C(3Fg + T{AT X + T*Bg- CBT^x > (13) 

where F is the known Bernstein operator for the system. 

For the non-dissipative cases, last expression reduces to the well-known Bernstein MHD 
energy principle 

^ S = \J l(*PPo£d 3 x + J {*(3Fg\d 3 x > (14) 
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from where the eigenmodes and eigenfrequencies are calculated as 

J = Jfrrtf* (15) 

and the stability criterion is obt ained by requiring the positivity of the potential energy of 
the perturbation (Galindo Trejo |l987| ) 

S 2 W p{Bernstem) = ^J C(3F^d 3 x (16) 

subject to the normalization condition that the total kinetic energy is equal to one. Thus, 
the dissipative principle and the new frequencies are respectively: 

S 2 W P = ^ J {£*(3Fg + T*AT t + T*B(- i*BT x )d 3 x > 0. (17) 



/ {CPF(+ T*AT X + T*B£-CBT*)d 3 



x 



^ = _ jvs . . ^i~s s (lg) 



with the same normalization condition. 



III. APPLICATION TO THE STABILITY OF A CORONAL INHOMOGENEOUS 
LOOP MODEL 

We are interested in analyzing the stability of non-homogeneous loops. This is, loops 
with inhomogeneous distributions of plasma density and temperatures. This character of 
the system poses additional difficulties for the stability analysis because the inhomogeneous 
nature of the medium determines the structure of the disturbance which is no longer sinu- 
soidal, making the traditional normal mode analysis useless for this treatment. Moreover, 
there may exist a continuous spectrum of stable modes besides the discrete one. As a first 
order approximation we neglect the effect of gravitational stratification and thus confine 
the analysis to characteristic spatial scales lower than the pressure scale height in the so- 
lar corona. In order to analyze the stability and to obtain the frequencies and modes the 
physical quantities in eq. [T7] and eq. [18] must be calculated along the loop structure. 



A. Mechanical equilibrium 



In order to determine an equilibrium configuration we assume force-free equations due 
to the fact that in plasma with low (3 (gas pressure over the magnetic pressure) the pressure 
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gradient can be neglected in comparison to the Lorentz force. The coronal arcade is obtain 
from the equations 



V x B = «B = (19) 

j x B = 0. (20) 

Also, Bo- Vp = and thus the pressure has a constant value along the loop. We assume that 
the unperturbed magnetic field is B = (B 0)X (x, z),0, B 0tZ (x, z)) and obtain the equilibrium 
field components 



B 0x = -Boocos(^-x)e ^ z (21) 
Bo* = fl M sin(^rz)e-&* (22) 



with 



The relation 



B = B 00 e-^ z e s . (23) 



2L. 
z = z t -\ in 

7T 



COSI7T— t 

v 2L' 



(24) 



is straightforward. The arc element s (see Figure [T]) can be expressed as 



ds = dx\j ( 1 + ^\ = dx A (25) 



with A = y/l + {z'f. 



B. Thermal equilibrium 

The thermal equilibrium is obtained from eq. [7] with L = (in eq. ED- Thus expressions 
F c = —koT^V \\T satisfies the two relations 

3F c dT —F c dF c p 2 TT . . 

oT os k T2 uT m 2 



12 



-L X =0 L 



FIG. 1: Schematic figure of the magnetic arcade with z(x) = zt + ^ln [cos(-7r^)], x and z the 
Cartesian coordinates, et and e n the tangential and normal versors respectively. T& and T t are the 
temperature values at the bottom and top respectively. The same notation is used for the density 
p. 



from where we obtain the two equations 



-dTkoT? 
ds 

— 2 Q{T') - H 



dT' 



(27) 



where we assume F c (s = 0) = as dT /ds = at s = and H = Hq, so the constant value 
of eq. His v = 0. We then can replace F c in eqs. E7| and give Q(T) its explicit expression. 
Then integrating between T t and Tf, (the temperatures at the top and the bottom of the 
loop respectively) and using (dT / ds)T=T b = and T t ^> T b we obtain the constant value of 
the heating function H = 7p 2 xT t a ~ 2 /(8k 2 B (a + §)). Also, we find 



dT 

ds 



P 2 X 



rpc 



T 



2-a 



(28) 



which is equivalent to the calculus in chapter 6 of Priest 1982J. Our aim is to obtain T 
as a function of the line element s. From eq. s = f(T) given as an integral expression 
of the temperature, which has to be inverted. Thus, for calculus purposes, we define v = 
1 — (T/T t ) 2 ~ a and replace T as a function of v in eq. [2SJ Then we obtain s = f(v) as 

1™ A 



where 



-1/ 



t) q - l dt 



(29) 



(30) 
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(Arfken and Weber Il995l ) with 



P 



1 

2' 



- + -)(2 

'2 4 A 



a) + 1, 



^ = (2 - a)T?~^ ((p 2 X )/(2k (a + jj)A&))*. 



Then, T = f 1 (s) as 



(is 



-4 



-i -l 



(31) 



C. The perturbation 

In order to calculate the stability and structure modes the general perturbation expression 
along the equilibrium loop is written 



6 = [C(x)e t + r](x)e n + £ y (x)ey]e 



iky 



(32) 



and 



7i = T^e^. 



(33) 



Then, representing the equilibrium functions of the different quantities with a sub-index 
and using the loop parameters and the mathematical relations presented in the Appendix, 
we obtain the non-dimensional expression for the energy principle (eq. [1] 



+(3T f 



^j^f + PoDxxf - k^p-iy - kp ^- 



+ 



A 



-z' n . 



14 



i . z' _ dB „ dri 



1 dx A y A diz; 



dx 

2 




+ 



t J- Z 

A +T] A 

dBn „ dri , z' ^ . 
^ + Bo ^- fc A B ^ 



+C 2 



i i rfT 0m dTi T n 2 d / 1 \ m dTi 

(^7^2 U^rn 1 9—— I I T -4- 

A 2 dx 1 dx A dx V A J 1 dx 



A 2 1 dx 2 dx V2"° A 2 dx 



+ 



dx 



+P Po (D x f-kQT 1 + 



1 dp ^ p a dTi z' don dTi , . m 



(34) 



A dx A dx A dx dx 

where the non-dimensional quantity S 2 W P / (xT® +1 p 2 L/m 2 ) replaces S 2 W P , and C\ = 

7 

BQ kBT t p t /mfi and C2 = c m 2 T t 2 / \L 2 T t ap 2 ) were used. From this variational principle 
we can then analyze stability and obtain the mode structure and the associated frequencies 
for the general mode given by eq. 



IV. RESULTS AND DISCUSSION 

In order to calculate modes and frequ encies we followed the schematic procedure de- 
scribed in Paper I and in Galindo Trejo 1987| . We used a symbolic manipulation pro- 



gram to integrate the equations. 5 2 W P and the perturbations were expanded in a three 
dimensional-Fourier basis that adjusts to border conditions. Thus, a quadratic form for 
5 2 W p was obtained and was minimized with the Ritz variational procedure. A matrix dis- 
crete eigenvalue problem subject to a normalization constraint was obtained. The procedure 
is equivalent to solving eq. [18] of our modified principle. Once the modes are obtained, the 
stability condition of eq. [17] for the generalized potential energy: S 2 W P > must be corrob- 
orated. The following values were used for the numerical calculation of the modes 



a = — - 

t 

(- -) 

l 2'5 J 



A 



_ 6 
~ 5 

5-*/ P 2 X 



-T, 3 ( 
2 * 



2k k 2 ^ 
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Coronal loop parameters: L = 10 10 cm (or L = lOOMm), = 10 4 K T t = 10 G K n e = 
10 8 cm~ 3 electron number density p t = 2kBT t ; pt = mp t /kBT t . 

Our main concern was to know whether the magnetic configuration of equilibrium could 
be stable under linear perturbations. For non homogeneous configurations it is well known 
that the stable eigenvalues ca n hav e continuo us sp ectra while the unstable ones have a dis- 



crete spectrum (see Freidberg [19821 ] or Priest 19821 ] ) . If the resulting mode components have 



a characteristic wavelength of the order of the equilibrium structure, the non-homogeneous 
character of modes could determine, for the stable modes, a continuous spectrum. Thus, in 
this case, the traditional normal mode analysis gives only a rough description because one 
of the consequence of the existence of the continuum is that an accumulation of discrete 
eigenvalues can take place at either boundary, generally at u 2 = or u 2 = oo, indicating 
the presence of a continuum stable spectrum. Note that as the basis used is discrete, the 
spectrum is necessarily discrete. However, the additional evaluation of the generalized po- 
tential energy provides the correct unstable modes and gives an approximate value of the 
most probable stable period when the smaller uj 2 is not located at the boundaries. 

We used different values for k: k = 0, k = 0.5 and k = 10 (k is the wavenumber associ- 
ated with the perturbation component transverse to the plane that contains the magnetic 
configuration). We also calculated the frequencies and modes for two different values of the 
magnetic field: B 00 — 11G and B 00 = 100C 

Table 1 and Table 2 show the eigenvalues (periods) associated with the different modes 
for the cases B 00 = 11G and -Boo = 100G respectively, considering k = and obtained by 
solving eq. [181 We obtained 12 eigenfrequencies and 12 eigenmodes for each of the magnetic 
field values i.e. we used a three-component expansion and a four-component perturbation 
vector. We evaluated the mode behaviour for k ^ 0. For each mode corresponding to a 
complex eigenvalue, the perturbation ^ y was at least two orders of magnitude smaller that 
the parallel ( component and the normal r] components. For the modes with real eigenvalue, 
in only one case was £ y comparable to the smaller of the two other spatial perturbations. 
Thus, for numerical simplicity, we used k = and we discarded three zero eigenvalues 
associated with this choice of k. We then analyzed only the 9 relevant modes. This means 
that a two-dimensional analysis of the dynamics of the problem is reasonably able to obtain 
the overall behaviour within the approximations we are considering. Thus, we decided to 
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FIG. 2: Mode components corresponding to the first mode Pi = 36.3 min for the cases: a) £ and 
B 00 = 11G] b) r) and B 00 = HG; c) ( and B 00 = 100G; d) r, and B 00 = 100G. 

investigate the unstable modes and to consider the most stable one as a reference value for 
stability. The most stable mode is the one that has a real uj value and minimizes 5 2 W P (it 
is the mode that gives the minimum positive value of the functional S 2 W P ) and the most 
unstable one is the mode that has a complex uj value with the minimum value of r ~ 1/M 
(r the instability time). 

From the analysis of the table data we can conclude: 1) for each of the two investigated 
magnetic values we have three complex values of uj and six real ones; 2) in the two magnetic 
field cases the eigenvalues of the first mode are the same; 3) in all the other cases the eigen- 
values with Boo — 100G are almost an order of magnitude smaller than the corresponding 
values of B o = 11G; 4) the series of eigenvalues is such that it could be possible that the 
stable periods accumulate at uj = 0, thus the definite stability characterization is subject to 
the evaluation of the generalized potential energy of the modes. 

We analyzed the structure of the modes with complex values of uj as they are possible 



candidates for instability (Freidberg |1982| ). We noted that in the two first modes the 
component that is tangent to the magnetic field | ( | is greater than the component | r\ \ 
that is normal to it. This can be seen from Figure [2] and Figure [3] where ( and r\ are shown 
for the cases: Boo — HG and Boo — 100G respectively, also using k = 0. The third mode 
(see Figure HI) has comparable values of | £ | and | r] \. 

The fact that for the first mode the two values of Boo give the same time eigenvalue 
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FIG. 3: Mode components corresponding to the second modes: a) £ component of P2 = 6.6 rain 
with £>oo = b) 77 component of P2 = 6.6 mire with Bqq = 11G; c) ( component of P2 = 0.7 min 
with i?oo = d) n component of P2 = 0.7 rain with £>oo = 100G. 





FIG. 4: £ and rj components for the third mode, a) left: P3 = 4.3 rain with Bqq = WG and b) 
right: P3 = 0.5 rain with Bqq = 100G. £ y has vanishing values. 

Pi = 2 ■ 60tt/uj = 36. 3i min indicates independence from the magnetic structure. This 
is consistent with the relative values between the two components in the two -Boo cases: 
I Q |3>| 77 I (see Figure [2]). Thus, these magnetoacoustic modes are more of the acoustic type 
I £ |^>| 77 I than of the Alfven type, i.e. | £ |<C| r\ \ (see Figure [5]). Also, the obtained period 
is included in a rang e (10mm < P < 60mm) where MHD slow acoustic modes are expected 



(Aschwanden [20041 ] )■ 

Figure [3] shows the second mode for B 00 = 11G and B 00 = 100G respectively. Also for 
both cases the | ( \ perturbation is greater than the normal perturbation | re | by an order 
of magnitude. 

Figure H] show the superposition of ( and re for the third modes corresponding to P = 
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FIG. 5: Schematic classification of fast and slow magnetoacoustic waves. 9 is the angle between 
the mode and the magnetic field: 9 = corresponds to large values of £ and 9 = tt/2 corresponds 
to large values of 77. 



4.3 min, B 00 = 11G and P = 0.5 min, B 00 = 100G respectively. Note that in these cases, 
when the component 77 is relevant, resembling an Alfven wave, the relation between the 
eigenvalues (periods) of the different magnetic fields is Pug — IOPloog, in accordance with 
the relation between the two values of B 00 (11G) ~ 10Poo(100G) and with the corresponding 
values of the Alfven velocities of the medium va = Poo/ yJJTp. 

Figure [5] gives a schematic classification of fast and slow magnetoacoustic waves from 
where we can analyze the behaviour of the modes. The first mode corresponds to 9 w and 
as its eigenvalue is independent of the magnetic field it gives a slow magnetoacoustic mode. 
The third mode corresponds t o < 6 < 7r/2 and as Pz,i\g — 10P3,ioog it looks like a fast 
magnetoacoustic mode (Priest [19821 ] ) . 

Then, in order to establish the final unstable modes we integrated eq. [T7| for each of the 
normal modes, i.e, the integrand is the generalized potential energy density. 

Figure [6] shows the generalized potential energy density as a function of the independent 
variable x for the three first modes (see Table 1), and for the most stable one which was 
P4. We show the case Poo = HG, the case with Poo = 100G has the same functional 
dependence. Table 3 shows the eigenvalues and the potential energy for the modes with 
complex eigenvalues and for the most stable one. Note that, even when u has complex values 
for the three first modes, as S 2 W P is positive in the second and third case, the Pi = 36.3 min 
mode is the only unstable one. The fact that, on the contrary to what happens with the first 
mode, the other modes with complex u seemed to accumulate at the origin is an indication 
of non-real unstable modes. 

Figure [7] shows the structure of the components ( and 77 for the most stable mode P4 
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FIG. 6: Generalized potential energy density as a function of x for the modes a) Pi = 36. min 
with 6 2 W P = -16.; b) P 2 = 6.6 min with 5 2 W P = 2.9; c) P 3 = 4.3 min with 5 2 W P = 6.63; d) 
P A = 3.4 min with 5 2 W P = 2143, (B 00 = 11G for all the cases). 

and for the two cases: B 00 = 11G and B 00 = 100G. Note that | rj \>\ C | and that 
Pa,b{ug) — 3-4 min ~ 10 • 0.36 min = 10 • P^b(ioog) 

The mode structure of the stable eigenvalues can also be compared with recent results 
from the literature. Magnetoacou stic os cillations of the fast kink type have been studied 



theoretically (Edwin and Roberts 1983J) and directly observed in EUV wavelength s with 



TRACE (an updated review of theoretical and observational results in Aschwanden 2004 1 
and references therein). The observations are usually modeled by cylinders with a surface 
boundary representing coronal loops. The dispersion relation is obtained matching the 
internal and external MHD solutions via the requirement of continuity of pressure and 
perpendicular velocity. As in our model, the observed kink-mode oscillations correspond to 
the long-wavelength regime. In coronal conditions the magnetic field is almost equal inside 
and outside of the loop and the kink oscillation speed is almost the Alfven one depending 
on the ratio of external and internal density values, i.e., outside and inside the loop. On 
the contrary, our model is performed by perturbing a magnetic arcade, without considering 
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FIG. 7: components of the most stable periods a) £ component of P4 
n component of P4 = 3.4 min with Bqq = 11G; c) £ component of P4 
d) n component of P4 = 0.36 min with Bqq = lOOG. 



= 3.4 min with Bqq = 11G; b) 
= 0.36 min with B 00 = 100G; 



a cylinder with different inside and outside conditions. In eleven observational kink-mode 
oscillations from w hich t he m agneti c field of the events can be inferred were obtained by 
Aschwanden et al. 2002| and 2003 ]. The comparison of our stable mode data Pi>i in the 
.Boo — case is in good agreement with the kink-mode observational results. The period 
range (see Table 1), the magnetic strength (_B o = HG) and the wave speed (Alfven speed) 
fit the observations for similar loop densities and loop lengths. Also, the stable modes Pj>i 
with £>oo = 100G (see Table 2) have periods that are comparable with the expected range of 
fast sausage- mode p eriods (P ± 1 — QOsec) and wave speeds of the order of the Alfven speed 
(Aschwanden 20041 ]). However -even when a more precise comparison requires a modeling 
that takes into account differences between external and internal conditions- it is worth 
investigating whether these type of modes could be associated with more intense magnetic 
fields in comparison to the associated kink-mode magnetic fields. This will be attempted in 
future work. 

A main result regarding stability is that the characteristic time r = 36 min in which the 
instability grows is large enough to guarantee a relative permanence of the structure before it 
fades away: r ~ t a f, s ; whe re t n b« is the typical characteristic time in which loops seem stable 
(see Costa and Stenborg 20041 ] ) . Thus, even when the non-linear stationary configuration 
of Figure [1] is unstable it lasts long enough for the observations to be made. Moreover, we 
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FIG. 8: Thermal perturbation (Ti component) for the cases: a) P = 36 min; Bqo = 11Gb) 
P = 36 min; B 00 = 100G 
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FIG. 9: Schematic description of the unstable mode superimposed on the magnetic structure. At 
a definite phase the perturbation is always positive, it grows until it reaches x = ±L/2, then 
decreases until it becomes zero at Z = Z t 

confirm that the dynamic brightenings usually observed could be due to magnetoacoustic 
waves i.e. the perturbations have short periods in comparison with the time that instability 
occurs: P 4 = 3.4 min and P 4 = 0.36 min satisfy P 4 r. 

Thus, even when further calculation is needed in order to adjust the characteristic times, 
it seems that wave-based models could be able to describe the scenario of non-isothermal 
coronal loops for sufficiently short times comparable with the characteristic time in which 
the instability grows and the structure fades away. A more speculative argument about the 
relation between wave-based models and flow-based ones is given in the conclusions. 
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FIG. 10: The curve formed by the resulting component perturbations in the vector space of per- 
turbations for a) Boo = HG and b) Boo = 100G. 

V. CONCLUSIONS 

We investigated -via a thermodynamic energy principle- the stability of a coronal inho- 
mogeneous loop model in a non-linear equilibrium state, i.e. a given thermal and magnetic 
equilibrium configuration. We also obtained the frequencies and their associated modes. The 
perturbation chosen was of the general type described by eq. [32] which allowed the study of a 
more complex mode structure with coupled thermal and mechanical displacements from the 
equilibrium state. We used a three-component Fourier basis expansion on the independent 
coordinate x to characterize the unstable modes. We obtained three complex eigenvalues 
and six real ones with their corresponding eigenvectors for each of the magnetic field values 
Boo analyzed. The other three modes were discarded. The definite stability condition of the 
modes is given by integrating the generalized potential energy density of eq. [TTJ, allowing the 
interpretation of long-wavelength disturbances that are present in inhomogeneous media. 

1. We used different values of A; (A; = 0, A; = 0.5 and k = 10) to calculate the eigenvectors 
with complex eigenvalues and for all cases we obtained vanishing values of t; y with re- 
spect to the other perturbed quantities. When we repeated the procedure to calculate 
the modes with real eigenvectors we obtained small but not vanishing values of £ y in 
comparison with the other components. Thus, two dimensional loop coronal models 
with a temperature gradient of two orders of magnitude between the bottom and top 
are a good approximation to study the whole three dimensional stability. 

2. We can classify the structure of the modes obtained as follows: a) those for which 
C 3> f] and b) those for which r\ > C- In the first case the perturbation V\ = d^/dt is 
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almost parallel to the magnetic field and the eigenvalue is relatively independent of its 
intensity resembling the acoustic waves where v s is independent of the magnetic field 
(see Figure E]). This basic longitudinal mode describes an oscillation between parallel 
plasma kinetic energy and plasma internal energy In the second case i>i = d^/dt 
has an important orthogonal component and the eigenvalue varies with the magnetic 
field P\\g — IOPioog resembling the dependence of the Alfven waves va — -Boo- When 
the wave is nearly transverse it describes an oscillation between perpendicular plasma 
kinetic energy and the combined magnetic compressional and line bending energies. 
Thus, the first case can be thought of as slow magnetoacoustic waves and the second 
one as fast magnetoacoustic waves. The period of the slow magnetoacoustic mode is 
also in accordance with observational data. Between the fast magnetoacoustic modes 
and in the long wavelength regime we distinguish two possible types, depending of the 
strength of the magnetic field. For the modes with Pi>\ and B 00 = 11G we found that 
the period range, the magnetic strength and the speed of the modes resemble a fast 
kink-mode. Also, through the consideration of period range values we suggest that 
modes with Pj>i and B 00 = 100G could be thought of as sausage modes. However, to 
go further with the classification of mode type a modeling that takes into consideration 
differences between the inside and outside of the loop is required. Also, the non- 
homogeneous character of the problem places serious limitations on conclusions in 
relation to stable modes. 

3. We found only one unstable mode with characteristic growing time: t u = 36 min. The 
approximate and most stable mode is P4 = 3.4 min for Poo = 11G and P4 = 0.36 min 
for B 00 = 100G. The fact that there is an unstable mode means that the equilibrium 
state is unstable and that wave-based models are not adequate to fit observations. 
However, as r u > P4 by an order of magnitude or two (depending on the Poo value) the 
equilibrium appearance of the loop and the brightening effects of the most stable mode 
could be sustained by a characteristic time which is in accordance with observational 
data (r n the characteristic time of the instability). 

4. A much more speculative argument, which needs further analysis and numerical cal- 
culation of the non-linear behaviour of the modes, is as follows: the non-linear growth 
of unstable modes influences the stable modes (they are called slaves), the resulting 
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behaviour is fundamentally governed by the most unstable modes. As we obtained a 
unique unstable mode, of the type of a slow magnetoacoustic wave, this indicates an 
overall unstable behaviour governed by the tangential ( component and the thermal 
one. The thermal component of the unstable mode is shown in Figure [SJ As a charac- 
teristic wavelength of the components ( and 7\ is L/2, it would be worth investigating 
whether this instability could be associated with a lim it-cycle solution genera lly char- 
acterized as a flow-based model (Gomez et al. 1990] ; De Groof et al. 2004] ). If this 
is the case, r should be the growth of the instability, before it reaches its non-linear 
saturated value, in a new equilibrium state of an oscillatory type in the ( and T\ 
components. Thus, both types of models (waves and flow) converge in explaining the 
instability of a magnetic structure with long wavelength perturbations of the order of 
the magnetic structure. Also, even when the modes were linearly unstable, the fact 
that the dominant varying components are T\ and (, with the last one parallel to 
the magnetic field could imply that the magnetic structure (but not the equilibrium 
state) lasts much longer than what is stated by r. Moreover, this is in accordance 
with the energetic description of the type of perturbation. Slow, nearly longitudinal 
magnetoacoustic modes describe a basic oscillation between parallel plasma kinetic 
energy and plasma internal energy where the magnetic energy plays no relevant role. 
This could justify long lasting loop observations with dynamic plasma inside. Figure [9] 
is a scheme of the unstable mode superimposed on the magnetic structure. In half of 
the period the perturbation is always positive and grows until it reaches x = ±L/2, 
then decreases until it becomes zero at x = and Z — Z t . The perturbation gives the 
tangential velocity V\ = d£/dt of the plasma particles at each point of the magnetic 
configuration. Thus, in half of the period, as described in the figure, the plasma is 
emerges from the chromosphere i.e., £ = (e n . In the other half, the perturbation is 
inverted with respect to the figure, it is always negative, i.e. £ = — C e n, and the plasma 
particles fall into the chromosphere. A limit cycle is known to be a closed curve (a 
cycle) in the vector space formed by the perturbations. Figure [HI] shows the result- 
ing curve in the space of perturbations (Ti,£) for the two magnetic fields studied. It 
seems that only for relatively large values of the magnetic field limit cycle are solutions 
possible. 
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VI. APPENDIX: MATHEMATICAL TOOLS 



The following equations and relations are needed in order to obtain eq. 



p ' = kM (35) 

dp dp Q ds dp dT dT 

A -7- -JZr = A -TT ( 36 ) 



dx ds dx ds dx ds 

with A = a/1 + (z') 2 . From Figure [1] it is easy to show 



and 



Gt ~ A + A Gz] e " A + A 6 ' 



Then, the spatial perturbation £ can be written in the Cartesian system as 

£= [f(C,v)e x + i£ y e y + g(C,v)e z ] J ky (37) 



taking into account 



(■ e x = ( ^ + ^ V (x) ) = /(C, v)e iky (38) 



and 



Then, 



and 



e, = ( ^ - ^(x) ) c *« = <?(C, ij)e*«. (39) 



<7(C,r7) = 77(^(1 + * -/(C^K) (40) 



„ „ d ( 1 \ . , . 1 rfC(^) d f z' ,z' d . , x 



26 



-Pi Pi P3 P4 P5 P& P7 Ps P9 Pw P11 P\i 



36. i 6.6 % 4.3 i 3.4 3.1 1.8 1.4 1.3 1.0 0.0 0.0 0.0 
TABLE I: Periods associated with the unstable and stable eigenvalues (minutes) for £>oo = HG. 

Pi Pi P3 Pa P5 Pq Pi P& P9 Pw Pu P\i 
36.3 i 0.7 % 0.5 i 0.36 0.33 0.2 0.15 0.14 0.1 0.0 0.0 0.0 



TABLE II: Periods associated with the unstable and stable eigenvalues (minutes) for B 00 = W0G. 

n „ d 2 ( 1 \ ^ , d 1 d((x) 1 d 2 . . 

rf 2 /V\ . , o d ( z'\ r)(x) z' d 2 , , 



Finally to obtain a non-dimensional equation the following changes were made 

P rj, T B x,z 

The other non-dimensional quantities are obtained immediately from these ones. 
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